//----------------------------------------------------------------------------
// M A R I T I M E  R E S E A R C H  I N S T I T U T E  N E T H E R L A N D S
//----------------------------------------------------------------------------
// Copyright (C) 2008 - MARIN - All rights reserved - http://www.marin.nl
//----------------------------------------------------------------------------
// Program    : mermaid
// Module     : vbm
// File       : FloatUtils.cpp
//----------------------------------------------------------------------------

#include <math.h>
#include <stdlib.h>
#include <algorithm> // min+max
#include "FloatUtils.h"
#include "assertplus.h"

#ifdef _WINDOWS
#pragma warning(disable: 4244) // conversion from 'double' to 'float', possible loss of data
#endif

static double ran2();

/// some necessary floating point routines
/** this namespace could be made into a singleton class when I have time,
    especially because of the initializations needed in the random generator */
namespace FloatUtils
{

/// raises base to exponent

/** (adapted from Jack Lyle) */
double power(double base, double exponent)
{
    const double cTiny = 1e-15;
    double result = 0.0;

    // Deal with most abundant cases
    if (exponent == 2)
    {
        return base * base;
    }
    else if (exponent == 3)
    {
        return base * base * base;
    }

    // Deal with the near zero special cases
    if (fabs(base) < cTiny)
    {
        base = 0.0;
    }
    if (fabs(exponent) < cTiny)
    {
        exponent = 0.0;
    }

    // Deal with the exactly zero cases
    if (fabs(base) < cTiny)
    {
        return 0.0;
    }
    if (fabs(exponent) < cTiny)
    {
        return 1.0;
    }

    // Cover everything else
    if (base < 0 && exponent < 0)
    {
        result = 1.0 / exp(-exponent * log(-base));
    }
    else if (base < 0 && exponent >= 0)
    {
        result = exp(exponent * log(-base));
    }
    else if (base > 0 && exponent < 0)
    {
        result = 1.0 / exp(-exponent * log(base));
    }
    else if (base > 0 && exponent >= 0)
    {
        result = exp(exponent * log(base));
    }
    else
    {
        ASSERTFALSE();
    }

    // Correct the sign
    if (base < 0 && fmod(exponent, 2.0) != 0.0)
    {
        result = -result;
    }
    return result;
}

/// modulo function as defined in fortran90 and as opposed to fortran 77's mod

/** c's fmod(x,y) behaves as fortran77's mod function which differs for negative arguments! \n\n
   MOD( 8, 5)  gives  3     MODULO( 8, 5)  gives  3 \n
   MOD(-8, 5)  gives -3     MODULO(-8, 5)  gives  2 \n
   MOD( 8,-5)  gives  3     MODULO( 8,-5)  gives -2 \n
   MOD(-8,-5)  gives -3     MODULO(-8,-5)  gives -3 \n
   From http://mathforum.org/library/drmath/view/52343.html : \n
   The mod function is defined as the amount by which a number exceeds the
   largest integer multiple of the divisor that is not greater than that number.*/
double modulo(double a, double b)
{
    int sign = 1;
    if (b < 0)
    {
        sign = -1;
        a = -a;
        b = -b;
    }
    double quotient = a / b;
    double integer;
    modf(quotient, &integer);
    double multiple = integer * b;
    if (multiple > a)
    {
        multiple -= b;
    }
    return sign * (a - multiple);
}

#define TwoPi 6.28318530717959

/// generates random number

/** gauss = true: according to Gaussian distribution (Delft Hydraulics GKl 940927) otherwise not */
double randomNumber(bool gauss)
{
    if (!gauss)
    {
        return ran2();
    }

    /* Initialized data */
    static bool set = false;
    static double N2 = 0;

    /* Local variables */
    double x, N1, U1, U2;
    double ampl, phase;

    if (!set)
    {
        set = true;
        U1 = ran2();
        //	while (U1 <= 0.0)
        //        {
        //	    U1 = ran2(); // komt nooit voor
        //	}
        U2 = ran2();
        ampl = sqrt(-2.0 * log(U1));
        phase = TwoPi * U2;
        N1 = ampl * cos(phase);
        N2 = ampl * sin(phase);
        x = N1;
    }
    else
    {
        set = false;
        x = N2;
    }

    return x;
} /* rnor1 */


} // namespace

/// Long period (> 2^1018) random number generator of L'Ecuyer
/// with Bays-Durham shuffle and added safeguards.
/// Returns a uniform random deviate between 0.0 and 1.0
/// (exclusive of the endpoint values).

#define NTAB 32
#define IQ1  53668
#define IQ2  52774
#define IA1  40014
#define IA2  40692
#define IR1  12211
#define IR2  3791
#define IM1  2147483563
#define IM2  2147483399
#define RNMX (1.0 - 1.2e-7) /// should approximate largest floating value less than 1

static double ran2()
{
    /* Initialized data */
    static int idum     = 0;
    static int idum2    = 123456789;
    static int iy       = 0;
    static int iv[NTAB] = { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 };

    /* Local variables */
    double ret_val;
    int j, k;

    if (idum == 0)
    {
        /* Initialize */
        idum  = 1;    /* Be sure to prevent idum = 0 */
    idum2 = idum;

        /* Load the shuffle table (after 8 warm-ups) */
    for (j = NTAB + 7; j >= 0; --j)
        {
        k = idum / IQ1;
        idum = IA1 * (idum - k * IQ1) - k * IR1;
        if (idum < 0)
            {
        idum += IM1;
        }
        if (j < NTAB)
            {
        iv[j] = idum;
        }
    }
    iy = iv[0];
    }

    /* Start here when not initializing. */

    /* Compute idum=mod(IA1*idum,IM1) without overflows by Schrage’s method */
    k = idum / IQ1;
    idum = IA1 * (idum - k * IQ1) - k * IR1;
    if (idum < 0)
    {
        //	idum += IM1; // never occurs (?)
    }

    /* Compute idum2=mod(IA2*idum2,IM2) without overflows by Schrage’s method */
    k = idum2 / IQ2;
    idum2 = IA2 * (idum2 - k * IQ2) - k * IR2;
    if (idum2 < 0)
    {
    idum2 += IM2;
    }
    j = iy / (1 + (IM1 - 1) / NTAB); /* Will be in the range 1:NTAB */

    /* Here idum is shuffled, idum and idum2 are combined to generate output */
    iy = iv[j] - idum2;
    iv[j] = idum;
    if (iy < 1)
    {
        iy += (IM1 - 1);
    }
    ret_val = std::min(iy * (1.0 / IM1), RNMX); /* Because users don’t expect endpoint values */

    return ret_val;
} /* ran2 */


